R Notebook

library(sf)
link='https://github.com/MAGALLANESJoseManuel/Midis_Bono600/raw/main/midisceplan.geojson'
midis=sf::read_sf(link)

names(midis)

 [1] "DEPARTAMENTO"                  "PROVINCIA"                    
 [3] "DISTRITO"                      "INSTITUCIO"                   
 [5] "UBIGEO_yana"                   "UBIGEO"                       
 [7] "programado_yana"               "ejecutado_yana"               
 [9] "yanaProp"                      "PeruLibre_congre"             
[11] "congreEmitidos"                "proporcion_congre"            
[13] "PeruLibre_presi"               "presiEmitidos"                
[15] "proporcion_presi"              "programado_b600"              
[17] "ejecutado_b600"                "nivel"                        
[19] "pobTotal"                      "pobDiscapacidad"              
[21] "superficie"                    "densidad"                     
[23] "altitud"                       "latitud"                      
[25] "longitud"                      "numeroCentrosPoblados"        
[27] "porcentajeDesnutCrónica5menos" "porcentajePobrezaExtrema"     
[29] "geometry"                     

str(midis)

sf [1,874 × 29] (S3: sf/tbl_df/tbl/data.frame)
 $ DEPARTAMENTO                 : chr [1:1874] "TACNA" "TACNA" "TACNA" "TACNA" ...
 $ PROVINCIA                    : chr [1:1874] "TACNA" "TACNA" "TACNA" "TACNA" ...
 $ DISTRITO                     : chr [1:1874] "CORONEL GREGORIO ALBARRACIN LANCHIPA" "POCOLLAY" "CALANA" "TACNA" ...
 $ INSTITUCIO                   : chr [1:1874] "IGN" "IGN" "IGN" "IGN" ...
 $ UBIGEO_yana                  : int [1:1874] 230110 230108 230103 230101 230109 230104 230303 180301 230102 230106 ...
 $ UBIGEO                       : int [1:1874] 230110 230108 230103 230101 230109 230104 230303 180301 230102 230106 ...
 $ programado_yana              : int [1:1874] 50230 7204 1879 40771 1766 21048 1076 26281 17153 1227 ...
 $ ejecutado_yana               : int [1:1874] 46451 6746 1753 37323 1625 19352 973 24819 15743 1133 ...
 $ yanaProp                     : num [1:1874] 0.925 0.936 0.933 0.915 0.92 ...
 $ PeruLibre_congre             : num [1:1874] 12877 1776 622 8565 353 ...
 $ congreEmitidos               : num [1:1874] 58524 10919 3117 72366 1664 ...
 $ proporcion_congre            : num [1:1874] 0.22 0.163 0.2 0.118 0.212 ...
 $ PeruLibre_presi              : num [1:1874] 19737 2845 943 12927 547 ...
 $ presiEmitidos                : num [1:1874] 59715 11168 3117 73341 1664 ...
 $ proporcion_presi             : num [1:1874] 0.331 0.255 0.303 0.176 0.329 ...
 $ programado_b600              : int [1:1874] 28770 4445 1218 26312 1240 11915 0 15989 10685 720 ...
 $ ejecutado_b600               : int [1:1874] 26662 4130 1158 24021 1171 10882 0 14975 9794 681 ...
 $ nivel                        : chr [1:1874] "Distrito" "Distrito" "Distrito" "Distrito" ...
 $ pobTotal                     : num [1:1874] 104353 17856 4441 130965 2747 ...
 $ pobDiscapacidad              : num [1:1874] 953 165 31 1436 9 ...
 $ superficie                   : num [1:1874] 188 266 108 1878 1116 ...
 $ densidad                     : num [1:1874] 555.84 67.22 40.98 69.74 2.46 ...
 $ altitud                      : num [1:1874] 562 690 881 583 397 ...
 $ latitud                      : num [1:1874] -18 -18 -17.9 -18 -17.9 ...
 $ longitud                     : num [1:1874] -70.3 -70.2 -70.2 -70.3 -70.6 ...
 $ numeroCentrosPoblados        : num [1:1874] 32 13 37 18 33 15 17 4 8 22 ...
 $ porcentajeDesnutCrónica5menos: num [1:1874] 3.3 2.81 3.33 2.76 3.73 ...
 $ porcentajePobrezaExtrema     : num [1:1874] 1.644 0.944 0.432 0.226 1.143 ...
 $ geometry                     :sfc_MULTIPOLYGON of length 1874; first list element: List of 1
  ..$ :List of 1
  .. ..$ : num [1:282, 1:2] -70.2 -70.2 -70.2 -70.2 -70.2 ...
  ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
  ..- attr(*, "names")= chr [1:28] "DEPARTAMENTO" "PROVINCIA" "DISTRITO" "INSTITUCIO" ...

h1=formula(ejecutado_yana~proporcion_presi+proporcion_congre)
reg1=glm(h1,offset = log(pobTotal),family = 'poisson',data=na.omit(midis))
summary(reg1)

Call:
glm(formula = h1, family = "poisson", data = na.omit(midis), 
    offset = log(pobTotal))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-220.060    -0.316     7.165    15.077    86.262  

Coefficients:
                    Estimate Std. Error  z value Pr(>|z|)    
(Intercept)       -1.2784567  0.0004798 -2664.37   <2e-16 ***
proporcion_presi   1.5801539  0.0075170   210.21   <2e-16 ***
proporcion_congre -0.1269199  0.0110913   -11.44   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1358372  on 1653  degrees of freedom
Residual deviance:  808249  on 1651  degrees of freedom
AIC: 824338

Number of Fisher Scoring iterations: 4

exp(coef(reg1))

             (Intercept)         proporcion_presi        proporcion_congre 
               0.9232442                1.1834246                0.8563463 
porcentajePobrezaExtrema 
               0.9994703 

library(ggplot2)
base=ggplot(data=midis,aes(x=b600prop,y=yanaProp)) 
base + geom_point()

Warning: Removed 899 rows containing missing values (geom_point).

sd(midis$ejecutado_yana)/mean(midis$ejecutado_yana)

[1] 2.531107

midis$propDisca=midis$pobDiscapacidad/midis$pobTotal

AER::dispersiontest(reg1,alternative='greater')$ p.value<0.05

[1] TRUE

h2=formula(ejecutado_yana~proporcion_presi+proporcion_congre+ porcentajeDesnutCrónica5menos+numeroCentrosPoblados+offset(log(pobres)))

na.omit(midis)

Simple feature collection with 888 features and 31 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -79.27475 ymin: -18.35067 xmax: -69.46805 ymax: -0.9609841
Geodetic CRS:  WGS 84

rbn=MASS::glm.nb(h2,data=na.omit(midis))#,control=glm.control(maxit=50))

Error in glm.fitter(x = X, y = Y, w = w, start = start, etastart = etastart,  : 
  NA/NaN/Inf in 'y'

rbn=MASS::glm.nb(h2,data=midis,control=glm.control(maxit=50))

summary(rbn)

Call:
MASS::glm.nb(formula = h2, data = midis, control = glm.control(maxit = 50), 
    init.theta = 4299.206932, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-6.4870  -0.4329   0.0480   0.4996   2.6237  

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -1.330e+00  3.429e-02 -38.790  < 2e-16 ***
proporcion_presi               2.621e-02  1.984e-02   1.321  0.18647    
proporcion_congre             -3.818e-02  2.714e-02  -1.407  0.15958    
porcentajePobrezaExtrema      -2.062e-04  8.739e-05  -2.360  0.01829 *  
porcentajeDesnutCrónica5menos -9.060e-04  1.059e-04  -8.558  < 2e-16 ***
numeroCentrosPoblados         -4.114e-05  1.540e-05  -2.672  0.00755 ** 
b600prop                       1.353e+00  3.707e-02  36.506  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(4299.207) family taken to be 1)

    Null deviance: 2296.33  on 919  degrees of freedom
Residual deviance:  733.12  on 913  degrees of freedom
  (954 observations deleted due to missingness)
AIC: 10191

Number of Fisher Scoring iterations: 1

              Theta:  4299 
          Std. Err.:  400 

 2 x log-likelihood:  -10175.07 

exp(coef(rbn))

                  (Intercept)              proporcion_presi 
                    0.2644901                     1.0265614 
            proporcion_congre      porcentajePobrezaExtrema 
                    0.9625415                     0.9997938 
porcentajeDesnutCrónica5menos         numeroCentrosPoblados 
                    0.9990944                     0.9999589 
                     b600prop 
                    3.8699414 

library(ggplot2)
dotwhisker::dwplot(rbn,exp=T)+ geom_vline(
           xintercept = 1,
           colour = "grey60",
           linetype = 2
       )

# + scale_y_discrete(labels=c("trabajo\nindependiente","analfabetismo")) + scale_color_discrete(name="Modelos para:\nCantidad de Asegurados")